Discrete stochastic models for traffic flow 



M. Schreckenberg^ , A. Schadschneider^ , K. NageP, and N. Ito^'* 

^ Institut fiir Theoretische Physik, Universitat zu Koln, D-50937 Koln, Germany, 
schreckOthp.uni-koeln.de and as@thp.uni-koeln.de 
^ Zentrum fiir Paralleles Rechnen ZPR, Universitat zu Koln, D-50923 Koln, Germany, 

kaiOzpr . uni-koeln . de 
^ CISC JAERI, Tokai, Ibaraki 319-11, Japan, ito@catalyst.tokai.jaeri.go.jp 

Abstract 

We investigate a probabilistic cellular automaton model which has been in- 
troduced recently. This model describes single-lane traffic flow on a ring and 
generalizes the asymmetric exclusion process models. We study the equilib- 
rium properties and calculate the so-called fundamental diagrams (flow vs. 
density) for parallel dynamics. This is done numerically by computer simu- 
lations of the model and by means of an improved mean-field approximation 
which takes into account short-range correlations. For cars with maximum 
velocity 1 the simplest non-trivial approximation gives the exact result. For 
higher velocities the analytical results, obtained by iterated application of the 
approximation scheme, are in excellent agreement with the numerical simula- 
tions. 
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113, Japan 
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I. INTRODUCTION 



Recently, there has been considerable interest in the investigation of traffic flow using 
methods of statistical physics. Independently in and cellular automaton models 
for the description of traffic flow have been proposed. Due to their computational simplicity, 
lattice gas automata have already successfully been applied to other problems, e.g. the 
simulation of fluids |Q (for further apphcations, see the book of Wolfram |^). 

A similar class of discrete models — which may be interpretated as traffic models or as 
models for surface roughening — have also been used for the description of the so-called 
asymmetric exclusion process (driven diffusion) [5-13]. Here several exact solutions have 
been obtained for processes where the particles can move at most one lattice spacing per 
update step. A natural generalization would allow particles to move over larger distances. 
This is more realistic with regard to modelling traffic since one usually has a whole spectrum 
of allowed car velocities. These generalized exclusion models are thus more appropriate for 



comparison with 'experiments' (i.e. measurements on freeway traffic fT^). 

In addition there have been a number of publications dealing with traffic flow in the 
framework of statistical mechanics, especially using cellular automata 15-17]. In most 
of these works two-dimensional traffic has been studied which corresponds to modelling the 
complex situation of traffic in the street-network of a city. Again the cars were allowed to 
move at most one lattice site per time step. It could be shown numerically [|^ that relaxing 
this constraint, i.e. the cars can have integer velocities up to an upper speed limit larger than 
one, behaves qualitatively as well as quantitatively in good agreement with real traffic (with 
an appropriate upper speed limit). It will be shown in this paper that, from an analytical 
point of view, the situation changes drastically when turning from maximum velocity one 
to higher speed limits. In this case real long-range correlations occur even in the stationary 
state which is not true for models with maximum velocity one. 

The exclusion models mentioned above may be classified according to the boundary 
conditions and dynamics used. There are two relevant types of boundary conditions: pe- 
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riodic and open. Using periodic boundary conditions one considers a ring on which the 
cars/particles can move ('Indianapohs situation'). Open boundaries correspond to the so- 
called 'bottleneck situation' where one imposes certain input and output flows at the chain 
ends. 

Basically one has to distinguish four types of dynamics. The dynamical variables may 
be updated one after the other in a certain order (sequential update), one after another in 
random order (random-sequential), in parallel for all sites of a given sublattice (sublattice 
update) or in parallel for all sites (parallel update). In the case of asymmetric update rules 
and sequential dynamics one has to distinguish at least two cases: update in or opposite 
to the direction of the traffic motion. Note that for open boundary conditions sequential 
and parallel dynamics are identical if the update direction is the same as the direction of 
traffic flow. For periodic boundary conditions, however, they are not identical. Here parallel 
dynamics correspond to sequential dynamics with a special site which closes the ring. This 
site memorizes a car even if it moved away one time-step before. This creates an obstacle 
in the ring giving rise to a lower flow through this site compared with sequential dynamics 
in the direction of the motion. 

In the present paper we use periodic boundary conditions and parallel update. In this 
case not much is known exactly, since all of the previous works [5-13] used either random- 
sequential or sublattice update. The advantage of parallel update with respect to sublattice 
or sequential update is that all sites are equivalent which should be the case in a realistic 
model with translational invariance. On the other hand parallel update enhances the possible 
system sizes in numerical simulations, especially because parallel or vector computers can 
be used easily. 

The paper is organized as follows: In Sect. 2 we introduce the model. In Sect. 3 the results 
of computer simulations are presented. We describe the different simulation techniques 
which have been used and compare their performance. In Sect. 4 we present a mean-field 
analysis of the model for arbitrary velocities Vmax- In Sect. 5 we introduce the so-called 
n— cluster approximation This improved mean-field method takes into account short- 
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range correlations. In the case Vmax = 1 we show that the 2— cluster approximation already 
gives the exact result. For Vmax = 2 we compare the results of the n— cluster approximation 
{n = 1, ... 6) with the computer simulations of Sect. 3 and find an excellent agreement. In 
App. A the exact solution for the case Vmax = 1 is derived whereas in App. B the equivalence 
of this case to an Ising model with repulsive interactions is shown. 

II. THE MODEL 

In the following we study single-lane traffic on a ring of length L with periodic boundary 
conditions. The model which has been introduced in is defined as follows: 
Each of the L sites can either be empty or occupied by one vehicle with velocity v = 
0, 1, . . . ,Vmax- At each discrete time-step t ^ t + 1 an arbitary arrangement of cars is 
updated according to the following rules: 

1) Acceleration: If the velocity w of a vehicle is lower than Vmax the speed is advanced 
by one [f = f -|- 1] . 

2) Slowing down (due to other cars): If the distance d to the next car ahead is not 
larger than v {d < v) the speed is reduced to d — 1 [v = d ~ 1]. 

3) Randomization: With probability p, the velocity of a vehicle (if greater than zero) 
is decreased by one [v = v — 1]. 

4) Car motion: Each vehicle is advanced v sites. 

These rules are applied to all cars in parallel (parallel update). The rules ensure that the 
total number A^ of cars is conserved under the dynamics (which is not true in the bottleneck- 
situation). Note that even for parallel update the randomization yields non-deterministic 
behaviour. For random-sequential update the probability p > is not essential because it 
only rescales the time axis |]T|. 
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In the simplest case Vmax = 1 the cars are allowed to move only one step during an 
update. For this situation several results are known [5-13]. Especially it can be shown 
that for random-sequential update the mean field Ansatz yields the exact equilibrium state 
which is equivalent to the fact that for a fixed number of cars every arrangement of 
cars occurs with the same probability. Therefore it is quite natural to take the mean-field 
approach also as a starting point for the investigation of higher velocities Vmax > 1 and 
parallel update. 

Our main interest will be the calculation of the so-called fundamental diagram (fiow 
q vs. density p = N/L). As described in [0] these results can be compared directly with 
measurements of real traffic [1^. One expects a transition from laminar fiow to start-stop 



waves with increasing car density. For Vmax = 1 it is easy to see that the fundamental 
diagram is symmetric with respect to p = 1/2 due to particle-hole symmetry. This is not 
true for realistic traffic where on finds a distinct asymmetry, i.e. the maximum of the fiow 
is shifted to lower values of p (~ 0.2). 

The rules given above take into account the basic features of real traffic. Firstly, they 
allow for a spectrum of velocities which seems to be necessary in order to break the particle 
hole symmetry of the model with only two velocities Vmax = 1- The maximum velocity 
still enters the system as a free parameter but it can be argued [|l| that most realistic are 
values around Vmax = 5. Secondly, the acceleration of the cars is very smooth compared to 
the deceleration which can occur in only one time step. For simplicity we assume the same 
maximum velocity for all cars but this condition is not essential. Finally the randomization 
step is necessary to avoid purely deterministic dynamics and to take into account natural 
fiuctuations of driving. 

The order of the update rules given above is crucial. Changing it would also change 
the performance of the model. The randomization step could e.g. be placed between accel- 
eration and deceleration (1-3-2-4) but then its infiuence is lowered drastically since every 
decelerating car does not feel any randomness. On the other hand one could also change 
the starting step. This will not infiuence the properties of the model but can simplify the 



calculations. If one begins with step 2 (as we will do for the n-cluster approximation) and 
proceeds 3-4-1 then one has the advantage that no cars with velocity zero occur since all 
cars were just accelerated by one unit. Therefore the number of possible states of a site is 
reduced by one. 



III. THE SIMULATIONS 

An integral part of our research were computer simulations. Simulations allow to obtain 
quantitative insight into a model in relatively short time. In this way, they complement 
the analytical work: Simulations are first used to test a variety of models until one has 
some overview and the most useful is found, and high quality simulation data are used 
to confirm analytical results (see Figs. |^ and |^). Further on, simulations are used to go 
beyond the analytically treated cases: either for variations/extensions of the model, or for 
more complicated quantities and issues such as the life-time distribution of the simulated 
traffic jams ||T^ or the possibility of self-organized criticality [Q. Last, but definitely no 
least, this work are the first steps towards an ultrafast microscopic simulation model for 
large traffic networks. We already have results for 2-lane-traffic ETl-pS and for real world 



network implementations including ramps, intersections, and junctions |pTip4|,p5 1 . 



A. Simulation technicalities 

Initially, we used a simple code on a workstation (and sometimes its replication on several 
nodes of a parallel computer) for getting an overview over the model's properties and for 
estimating its relevance for real world traffic [|l[]. Later, we implemented vectorizing and/or 
parallelized versions in order to obtain faster simulation speeds. 

Besides the advantage of getting high quality data in relatively short time, we have 
taken advantage of the simplicity of the model to implement it with different algorithms on 
many different supercomputers. This gives us intuition on how to implement more complex 



models p6[, and to make predictions on how fast our simulations will be in comparison to 



other microscopic traffic simulation models p^- the cases of 2-lane-traffic and of the 



network implementation, our predictions were quite accurate. 

For the practical coding, we considered three different approaches: site oriented, particle 
oriented, and an intermediate scheme. Site oriented directly implements the CA: A street is 
represented by a chain v of integers with values between —1 and Vmax, where —1 means that 
there is no particle at this site, whereas the other values denote a particle and its velocity. In 
contrast, particle oriented means that two lists {xi)i=i^,„^N and (fi)i=i,...,Ar contain position 
Xi and velocity Vi of each particle i {i = 1, . . . , N). This is similar to a molecular dynamics 
algorithm, except that particles are constrained to integer positions and velocities. 

Obviously, the particle-oriented approach will always be faster than the site-oriented 
one for sufficiently low vehicle densities. The particle-oriented approach is more flexible, 
and since in single-lane simulations passing of vehicles is not possible, the particle lists are 
always ordered, making efficient codes for all kinds of computers easy to write. In addition, 
an extension to continuous position and velocity is straightforward p8| . 

On the other hand, for the site-oriented (cellular automaton) approach single-bit cod- 
ing P3 is possible. This means that the model is formulated in logical variables, which 



may be stored bitwise into computer words. Logical operations on computer words treat all 
bits of the word simultaneously, giving a theoretical speedup of b, where b is the number 
of bits per word (usually 32 or 64). However, the practical gain for traffic simulations on a 
workstation is much lower because the bit-oriented approach cannot take advantage of the 
fact that only a fraction of all sites is occupied by a particle. Nevertheless, we found that, 
on a workstation, the single-bit algorithm is faster for densities above 0.05 (for Vmax = 5). 
In addition, the single-bit code runs very efficiently on a Thinking Machines CM-5 using 
dataparallel CM- Fortran and on a NEC-SX/3 traditional vector computer. The simulation 
data for the fundamental diagrams have been obtained this way. 

Once passing of vehicles is allowed (multi-lane traffic), for the particle oriented approach 
efficient memory allocation for parallel and/or vector processors becomes more difficult, and 
single-bit coding for the site oriented approach becomes tiresome. These observations led 



to a third, intermediate approach. As in the site oriented approach, each site is in one of 
{vmax + 2) states, but for the update only the relevant sites are considered. It turns out (see 
below) that on parallel but not vectorizing computers this algorithm is about as fast as the 
single-bit version. 

In Table |, we give an overview of the computational speeds on selected computers (see 



27| for more details on most of these results). All values are valid for a vehicle density of 
p = 0.1 and system size of L = 1333 333 sites, corresponding to 10 000 single lane kilome- 
ters. MUPS corresponds to Mega- Updates Per Second, i.e., the number of sites updated 
per second divided by 10^. These values are useful in order to compare with other imple- 
mentations of similar cellular automata or particle hopping systems. The other number is 
the extrapolated real time system size, i.e. the extrapolated system size (assuming linear 
speed-up) where the computation would be just as fast as reality. 

The most noteworthy features of the table are emphasized in boldface: 

• On vectorizing computers such as the SX-3 and the CM-5, the single-bit algorithm has 
a notable advantage over the other algorithms. On all other machines, the intermediate 
algorithm is never more than a factor 2.5 slower than the single-bit algorithm. 

• Already on a relatively modest machine such as an Intel Paragon with 64 nodes, our 
real time limit is 280 000 single lane kilometers. Since, e.g., the whole freeway net- 
work of Germany is about 60 000 single lane kilometers long (12 000 km x 2 directions 
X 2.5 lanes), the real use of this computational speed will be (i) real time applica- 
tions, where the traffic forecast has to be computed before the situation arrives; and 
(ii) Monte Carlo simulations, where many runs are necessary. 

• The by far fastest "realistic" traffic microsimulation world-wide to date is the PARAM- 
ICS microsimulation project |^. Their real time limit is ~ 20 000 km on 16k CM-200, 
a machine which is in terms of peak performance slightly faster than a 128-node- 
Paragon. In other words, it seems that our not completely realistic car following logic 
buys us about an order of magnitude in performance. 
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B. Simulation results 



Figs. |T| and |^ show typical time evolutions of different versions of the model. Fig. ^ shows 
parallel update with Vmax = 1 on the left, and random sequential update with Vmax = 5 on 
the right, both at density c = 0.5. Random sequential update for Vmax = 1 looks the same as 
the left picture. Obviously, neither using parallel update for Vmax = 1 nor taking a Vmax > 1 
for random sequential update does change the phenomenlogical behavior from the stochastic 
asymmetric exclusion process. For Vmax = 1, P = 0.5 corresponds to maximum flow, and 
in consequence, wave structures do not move in space For Vmax > 1, the point of 

maximum throughput is shifted to lower densities, and in consequence, in the left plot the 
waves are moving backwards. 

The situation is different when one combines parallel update and Vmax > 1 (Fig- Here, 
in the regime of maximum throughput (left plot) waves are only sparse, and they clearly 
move backwards. And even at higher densities (right plot), waves are much more distinct 
than in the random sequential case. 

In short, one can divide the models between random sequential update with arbitrary 
maximum velocity on one hand, and parallel update with maximum velocity Vmax > 2 on 
the other hand. Parallel update with Vmax = 1 phenomenologically is an intermediate case. 
Interestingly, it will turn out that this structure is reflected in the analytical calculations 
below: For random sequential update, the mean field solution is already exact. For parallel 
update and Vmax = 1, the situation is only slightly different because already the first step 
beyond mean field is exact. The situation is completely different for higher velocities in 
connection with parallel update, where the analytical approximations only converge slowly 
towards the simulation result. 

In addition. Fig. ^ allows an interesting comparison with fluid-dynamical models. Start- 
ing from ordered initial conditions, one clearly sees how instabilities develop and produce 
the start-stop- waves, very similar to results in Working out these connections is the 
topic of current research [p5| . 
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Fig. ^ (a) shows current vs. density curves for maximum velocities Vmax between 1 and 
5, plus for a different fluctuation parameter p = 0.25 at v^ax = 5. One clearly sees that 
the maximum throughput increases with increasing Vmax, whereas the density of maximum 
throughput decreases. In reality, the density of maximum throughput lies between p = 
0.15 and p = 0.2; it is given by the maximum speed of trucks which dominate the speed 
distribution for traffic near capacity [^. For that reason, having a higher speed limit for 



passenger cars does not help throughput; in many cases, it actually makes things worse [[33 

However, reducing the fluctuation parameter p increases throughput enormously. This 
is mostly due to the better acceleration behavior in that case [Q. 

In general, one sees that by varying the parameters f max and p, the fundamental diagram 
can easily be adapted to real traffic situations, although some of the underlying vehicle 
dynamics remain somewhat unrealistic: E.g. average acceleration from to 100 km/h takes 
place in 10 seconds. That fact that the fundamental diagram is nevertheless quite realistic 



is due to the fact that the first time steps fo the acceleration matter most And here, 
4 seconds for an acceleration from to 40 km/h are far more realistic. 

Another quantity of interest for traffic engineers are the velocity fluctuations 



0-[Vloc) ■= V {'^loc - VlocY , 



where vioc is the "local" velocity of vehicles crossing a fixed line. (The average of the lo- 
cal velocity is different from the usual ensemble average. For example, cars with velocity 
zero never enter the local average.) According to measurements and fluiddynamical argu- 



ments [^], these fluctuations are a good indicator of traffic near capacity. And indeed do we 



find for our model (Fig. ^ (b)) that these fluctuations very abruptly increase near capacity. 

IV. MEAN-FIELD THEORY 

The complete analytic solution of the traffic model is not possible except for the case of 
maximum velocity fmax = 1 (see Sect. 5.1 and Appendix A). Therefore some approximation 
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is necessary when one tries to study this model analytically for Vmax > 1- The following 
mean-field type theory will be the first step of such an analytical approach. 

The calculation starts from the stochastic description of the traffic model. Instead of 
specifing the car position and their velocities, we analyze the probability distribution of each 
site at each time step. We denote the probabihty that there is no car at site i (i = 1, 2, 3, 
• • •, L) at timestep t by d{i, t) and the probability that there is a car with velocity a {a = 0, 
1, 2, • • •, Vmax) at site i and timestep t by Ca{i,t). The c and d distributions together take 
into account all possible states of the different sites. Therefore one has the normalization 
condition for all sites and all timesteps 

d{i, t) + Co(i, t) + ci(i, t) + C2{i, t) + C3(i, t) + ■ ■ ■ c,_(i, t) = l. (4.1) 

Denoting with c{i,t) the total probability for site i to be occupied at timestep t, i.e. 
X]j=cr Cj{i, t), one simply has d{i, t) + c(i, t) = 1. 

According to the update rules in four stages (see Sect. 2) the time evolution of these 
probability distributions can be described by the following four sets of equations: 

• Acceleration Stage 

Co{i,ti) = 

Caii, h) = Ca-i{i, t) , <a < V^a^ (4.2) 



• Deceleration Stage 

co{i, h) = co{i, h) + c{i + 1, h) ^ Ci3{i, h) 

13=1 

a Umax 

Ca{ht2) = c{i + a + l,ti)X{d{i+jM) E c/3(i,ti) (4.3) 

j=l /3=a+l 

a 

+ Ca{i,ti)J\d{i+ j,ti), < Q; < Vmax 
3=1 

''^max 
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• Randomization Stage 

Co{i,t3) ^ Co{i,t2) +pci{i,t2) 

Ca{i,t3) = qCa{i,t2) + pCa+l{i,t2), < Qi < Wmax (4.4) 

• Motion Stage 

Ca{i,t + 1) ^ Ca{i - a,t3), < Q; < f max (4.5) 

The time t is assumed to take on only integer values, ti, t2 and denote the intermediate 
time steps between t and t+1 (sometimes defined asi+1/4, i + 1/2 and i + 3/4, respec- 
tively) . These time-evolution equations conserve for periodic boundary conditions the total 
number of cars at each stage . This formulation of the dynamics neglects spatial correlations 
completely since one assumes that all expectation values factorize into local terms. 

The variables c and d are real valued between and 1. The stochastic description 
originates from the stochastic nature of the third randomization step. The other three steps 
are purely deterministic. 

These time evolution equations are non-linear and further analysis of the full system has 
not been successful up to now. However, in the stationary state, i.e. in the limit t — > cxd, the 
c and d distributions become homogeneous in space (for periodic boundary conditions) and 
therefore, apart from the time dependence, also the site dependence can be omitted. Using 
this and combining the four update steps one gets the following set of v^aax + 1 equations: 

Co = (c + pd)co + (1 + pd)c Cf3 

13=1 

Ca = d°'[<lCa-i + {qc + pd)ca + {q + pd)c J2 (4.6) 

< a < t^max - 1 

Cwx-1 = d'''^"~^[qcv^..-2 + {qc + PC?)(c„^a.-l + C^^J] 

Climax qd ^ [Cllmax — 1 "I" Cumax] 
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These equations essentially describe the motion of a single car with statistical (densitity 
dependent) 'obstacles'. A remarkable feature of the equations is that they are linear when 
one specifies the density c = 1 — of cars. Therefore (4.6) can be recast in matrix form 
as Ac = c. The matrix A can be read off from (4.6), c is the vector with elements Cq,, 
a = 0, ...,Vmax- For small v^ax we can invert A to find the densities Cq, explicitly. For large 
values of Vmax it is more convenient to write down a recursion relation in order to obtain 
the steady state solution. 

From the first equation in (4.6) one can determine cq directly without knowledge of the 
other c'^s to give 

^0 = c ^3^. (4.7) 

Using this result and the second equation of (4.6) one can also write down the expression 
for ci 

2 1 + d + pd'^ 

For a larger than one a recursion relation can be derived calculating Cq, — dca-i using again 
the second equations of (4.6) 

a 



l + (g-p)ci" ^ qd' 



for a = 2, 3, ■ ■ ■, Wmax — 2. Therefore, starting with the expressions (4.7), (4.8) for cq and Ci, 
one can estimate C2, C3, ■ ■ ■, recursively. These results do not depend on the actual 

value of Vmax and thus are valid generally (provided Vmax > 2). 
Finally, from the last two equations of (4.6) one gets 



LTiiax 



qd'' 

•^■ymax ~ ~^ '17, '^fmax-1- (4--'--'-) 

1 — qa"^'^^ 

The Wmax dependence only occurs in these two quantities. In Fig. ^ some of the densities are 
shown for large fmax- The densities of the fast cars go to zero rapidly, since one expects an 
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exponentially fast decay from the recursion relations, (4.9)-(4.11). The contributions from 
cars with high velocities therefore are neglegible. We will mainly be interested in the flow 
f{c,p), defined by 

/(c,p)= 5:«c,. (4.12) 

a=l 

In the limit of f max oo it is possible to carry on the analysis using the iteration equations 
(4.7)-(4.9) and the generating function 

oo 

giz) = Y.z''co^. (4.13) 

As usual one simply has g'{l) = f{c,p). The equations (4.7)-(4.9) now can be combined to 
give one single equation for g{z) 

g{z)[l-dz] - g{dz)(f{l- z)[p + qz] = + pc^d{l - z). (4.14) 

Since fmax is infinite one does not have to worry about the upper boundary equations (4.10) 
and (4.11) for 

Ct)max-i '^vmnLx whosc coutributlon is neglegible. One should notice that 
the generating function g occurs with two different arguments {z and dz) which makes it 
impossible to solve this (linear) equation for g{z) explicitly. 
After differentiation an expression for g'{l) is obtained 

g'{l) = d{l~pc)-g{d)- (4.15) 

c 

The problem therefore is reduced to find an expression for g{d). This can be found by 
successive application of equation (4.14) with z = d^, n = 1, 2, 3, . . .. The final result is an 
asymptotic expression for g'{l) = f{c,p) 



f{c,p) = qcd 



n-l 



n=l 1=0 



(4.16) 



In Fig. ^ the fiow / is shown as a function of the concentration c of cars for various deceler- 
ation probabilities p. The slope at the origin (c ~ 0) of the fundamental diagram is infinite 
whereas the slope for c ~ 1 is —q. This mean-field result yields, compared with the simu- 
lation data shown above, much too small values of the fiow. This can easily be understood 
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since the reduction to a single car problem ignores all spatial correlations of the cars. Cars 
with high velocities tend to be equidistant and can therefore maintain a high velocity with a 
larger probability than in the mean-field system where it is much more difficult to accelerate 
and stay at high velocities over a certain period. 

Furthermore, even for Vmax = 1, the mean- field solution does not coincide with the exact 
result shown below. Note that for random sequential update the mean-field solution is exact 

for Vmax = 1- 



V. IMPROVED MEAN-FIELD THEORIES 

In order to improve the simple mean-field theory of the preceeding section we have to 
take into account correlations between neighbouring sites ||18|| . 

We divide the lattice into segments or clusters of length n {n = 1,2, . . .) such that two 
neighbouring clusters have n — 1 sites in common. The probability to find in equilibrium 
a cluster in state (ai, . . . , (T„) will be denoted by P„(cri, . . . , o"„). Due to the translational 
invariance of equilibrium state of the system with periodic boundary conditions one has not 
to specify the actual location of the n-spin cluster. In order to simplify the calculations 
we apply, as mentioned above, the four update-rules in the order 2-3-4-1 instead of 1-2- 
3-4 |jl8|. This has the advantage that after one update cycle one ends up with step 1 and 
therefore no car has velocity 0. It follows that every site j is in one of the Vmax states 
(jj = 0, 1, . . . , Vmax where now denotes an empty site. This change in the ordering finally 
has to be taken into account in the calculation of the fiow. 

The equilibrium probabilities for an n-site cluster are then determined by 

{rj} 

where the probability P2v^ax+n{r-vmax+i^ ■ ■ ^n+vmax) for clusters of length 2vmax + n has to 
be expressed through the n-cluster probabilities P„(ri, r„). This enlargement of the clus- 
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ter length occurs since all cars which can drive into or out of the cluster (cxi, . . . , (T„) within 
the next timestep contribute to the transition rates W^({o"j}|{Tj}). Thus we have to take 
into account not only the given cluster but also the Vmax sites to its left (with state variables 
'T-v^a^+i, ■ ■ ■ .tq) and the Vmax sites to its right (with state variables r„+i, . . . , r„+^„^J. The 
decomposition of the {2vmax + n)-cluster can be carried out by introducing the conditional 
probabilities 

P[Ti\ Ti+i,...,Ti+n-l ) = ^ p . T (5.2) 

at the left side and 

r)/ \ \ Pn{Tii "^i+li ■ ■ ■ i Tj+n-l) q\ 

r^[ Ti, ■ ■ ■ , Ti+n-2 \Ti+n-l) — — r [0.6) 

at the right side. With this definition we rewrite P2v^^^+n (in the n-cluster approximation) 
in the following form: 



P2Vmax+n{,T_Vr„ax + ly • • • 5 '^n+Vmax) = II (^i | i ••■^ 'Tj + n-l ) X (5.4) 

max 

X Pni'Tl, • • • ; Tn) H P{ 'Ti+l, Ti+n-l lTj+n) ■ 
i=l 

The transition probability W{ai, . . . , cr„|r_^^^^+i, . . . , Tn+v„^ax) depends on the probability p 
and vanishes if the configuration {r^vmax+iy • ■ ■ 5 '^n+v,nax) cannot evolve in one timestep into 
((Ti, . . . ,cr„) according to the rules 1) — 4) of Sect. 2. If W is non-zero it is of the form 
pniqn2 lutegcrs ni,n2 describing how many cars have to be decelerated (rii) through 

the randomization step when total number of cars which can drive is rii + n2. With the 
approximation ( |5.4|) it is possible to write down a closed system of equations for the n- 
cluster probabilities P„(o"i, . . . , cr„). The number of the equations is given by {vmax + 1)", 
the total number of possible configurations of n site variable with Vmax + 1 possible statesQ 
(without change of the order of the update steps one would have {vmax + 2)" equations). 



^In practice some of these equations turn out to be trivial so that the relevant number is less than 

{Vmax ~l~ 1) ■ 
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Due to the exponential growth with respect to n one is, especially for larger Vmax, restricted 
to only small cluster lengths n (For the realistic value of v^ax = 5 one has for the 2-cluster 
approximation already 36 equations!) 

The above approximation converges for n ^ oo to the exact solution for an infinite 
system (i.e. in the thermodynamic limit L oo. This approximation scheme is well known 



in the literature. It is analogous to the {n,n — 1)— cluster approximation of pBf, the n — 1 



step Markovian approximation of or the local structure theory of It goes back to 



the probability path method introduced by Kikuchi . 

With the knowledge of the n-cluster probabilities Pn(cri, • • • , o'n) it is then easy to calcu- 
late the fundamental diagram, i.e. the flow / as a function of the density c of cars. Since we 
have changed the order of the update steps one has to take this into account by performing 
the steps 2-3-4 at the end since the last step must be 4 (= car motion). After that one 
simply can calculate the density Ca of cars which will drive a sites in the next timestep 

Ca= Yl Pn{a,(^2, ■ ■ ■ ,(7n) (5.5) 

CT2,...,crn 

and then one proceeds as in the mean- field approximation eqn. (4.12) and calculates / = 



A. Vfjiax — 1 

In the case Vmax = 1 the site variables take on only the values a = 0, 1 where a = 
means no car and a = 1 a car with velocity one. In the 2-cluster approximation (( |5.1|) with 
n = 2) one has, according to the above arguments, a system of 4 equations. This can be be 
reduced to only one equation for P(l, 0) very easily with the help of the relations 

P(1,0) = P(0,1) 

P(0,0) = l-c-P(l,0), (5.6) 
P(l,l) = c-P(l,0). 
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The first equation is due to the 'particle-hole' symmetry P(l, 0) = P(0, 1) (in a closed ring 
one must have the same number of (0,1) and (1,0) pairs, therefore occuring with the same 
probability). The other two equations describe the conservation of cars in the system. The 
remaining equation for P(l, 0) reads 

gp2(l,0)-P(l,0) + c(l-c) = 0. (5.7) 

In the thermodynamic limit we therefore find [IS 



1 - Jl - 4gc(l - c) 

P(0, 1) = P(l, 0) = ^ ' ' . (5.8) 

Zq 

Going to the 3- and higher-cluster approximations one finds that the solution remains the 
same indicating that this is the exact result. In App. A we indeed prove that the solution 
( p.8| ) is exact in the thermodynamic limit. For finite systems the prove is also valid but one 
has to take into account an additional correction (normalization) due to the constraint of a 
fixed number of cars. 

The correct result, valid for any system size, is is given by 

P(Ar,L) = l5:'nP(a,,a,+i) (5.9) 

where M denotes the normalization constant and the sum Y^' runs over all configurations 
with a fixed number A^ of cars (i.e. Y^^=i = A^)- 

In contrast to random-sequential dynamics parallel dynamics leads to an effective at- 
traction between 'particles' and 'holes' (i.e. P(0)P(1) = c(l — c) < P(0, 1)) and thus to 
a higher flow. In App. B the mapping of this model to an equivalent Ising-model with 
antiferromagnetic next-nearest-neighbour interactions and an nonvanishing external field is 
shown. Due to the antiferromagnetic interactions the system shows an effective attraction 
of cars over two lattice sites, thus taking into account the well known effect of car bunching 
or platooning [pO|-|i2]. 



In order to calculate the flow one first has to perform the steps 2-3-4 yielding the 
probabilities P(t) to find a site in the state r = —1,0,1 (where now an empty site is 
denoted by r = — 1 ) after the fourth step of the updating procedure: 



P(-l) = l-c, P(0) = gP(l,l), P(l) = gP(0,l) (5.10) 
Therefore the fow is finally given by 



I - Jl-Aqcil-c) 
f{p, c) = qP{l, 0) = V ^_ ^5 

In Fig. 6 the exact result for the flow is shown (from the 2-cluster approximation) in com- 
parison to numerical simulations and the mean- field approximation {p = 1/2). One can see 
that the numerical and analytical data are in excellent agreement. In the deterministic case 
p = the flow is a linear function of c: / = (1 — |1 — 2c|)/2 |]T8| . 
The mean velocity per car v is then 

v = -Y.^nr) = P{l)/p. (5.12) 

P r=0 

In the deterministic case p = 1 this simplifies to 

1 for0<c<l/2 

(5.13) 

(l-c)/c forl/2<c<l 
which is the result of [01. 



The case Vmax = 2 is qualitatively very different from the case v^ax = 1- The flow 
diagram is no longer symmetric around c = 1/2. The n— cluster approximation seems not 
to become exact for any finite n, i.e. in this case long-ranged correlations are important. 
We have calculated the fundamental diagram for n = 1, . . . , 5. As shown in Fig. 7 the 
approximation converges fast to the results obtained by simulations. The difference between 
n = 4 and = 5 is less then 1%. 

The observation that the approximation does not become exact for small n reflects the 
fact that the physics for Vmax > 2 is distinctly different from the case Vmax = 1. As explained 
in Sec. 2 the regime looks qualitatively different both from Vmax = 1 (arbitrary update) as 
well as from random sequential update (arbitrary Vmax)- Moreover, in the case Vmax > 2 
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(parallel update), jams show characteristic branching behavior which is not observed for 

Vmax = 10- 

VI. CONCLUSIONS 

We have introduced and investigated a statistical model capable to describe accurately 
real traffic. Through the introduction of higher velocities it was possible to produce the 
asymmetric fundamental diagrams and the characteristic start-stop waves observed in real 
traffic. Through simulation, different regimes depending on the type of update and the max- 
imum velocity have been identified. Realistic for traffic is a model with parallel stochastic 
update and a suitable choice of the maximum velocity Vmax = 5. Already for Vmax = 1 
the stationary state is more complicated than for random sequential update. An effective 
'antiferromagnetic' interaction between cars favors equal spacing and in consequence higher 
throughput. 

Taking into account two-point correlations the case v^ax = 1 can be solved exactly. 
This is no longer true for higher maximum velocties where correlations become long ranged. 
Nevertheless, it could be shown that the n-cluster approximation converges fast to the 
simulation data. 
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APPENDIX A: EXACT SOLUTION FOR Vmax = 1 



In this Appendix we show that the stationary state of the model with Vmax = 1 is in fact 
given by equation (5.6) with an appropriate normahzation constant Z. The complete set of 
evolution equations for parallel update reads: 

Pt+iii^}) = E W{{ct} I {r}) ■ P,({t}), (A.l) 

{T} 

where Pt({cr}) denotes the probability for state {cr} = {o"i, . . . , a^} at time t. The transition 
probability W^({cr} | {r}) to move in one timestep from state {r} to state {cr} factorizes 



into local terms: 



L 

W{{(t} I {r}) = n Wia,, | r„ r,+i). (A.2) 



i=l 

The only non-vanishing transition probabilities are given by: 



W^(0,1|1,0) 


= 1 — p 


iy(i,o|i,o) 


= p 


W{l,a'\a,l) 


= 1 


iy((j,o|o,o) 


= 1 


W{0,a\0,l) 


= 1, 0": 



(A.3) 



whereas W{t, t' \ a, a') = in all other cases. (Note that in eqn. (5) of ref. [|l| due to a 
misprint a factor (1 — ai^i)/2 is missing on the right-hand side.) We now make the Ansatz 
that the probability in the stationary state P({cr}) factorizes into local two-site terms Po-io-i+i : 

P{{^}) = i{P...... (A.4) 

with periodic boundary condition (Tl+i = cti- We will define no-^a-'(cr) as the number of pairs 
of next neighbours (a, a') in a particular state cr. Due to the particle-hole symmetry of the 
system one has noi = niQ and the following simple relations for the system size L and the 
(conserved) number of particles hold: 
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L = 2 • noi + nil + rioo (A.5) 
N = noi + nu. (A.6) 

In the stationary state equation (A.l) becomes time-independent and the states {r} in the 
summation on the right-hand side can be classified according to the number of particles I 
which have to be moved in order to obtain the new state cr: 

^^5(noi,/,A)(l-p)y°i-'-^(PoiPio)"'^^~^Pir+^n^o"'"°'~""+^- 

1=0 A 

The summation index A is defined as A = nii(T) — nii(cr). Therefore the range of the A- 
summation depends on the particular state <t. The function g{noi, I, A) counts the number 
of possible states r leading to state cr given fixed values of rioi , / and A and reflects a kind 
of degeneracy. Summing g over A yields simply 

E^Ki,/,A)=(^°^) (A.8) 

Dividing eqn. (A. 7) by the left-hand side one gets 

noi / P P \^ 

1=0 A \pninoj 

Setting the expression in the parenthesis on the right-hand side to 1, (A.9) reduces with the 
help of (A.8) to an identity. Therefore the condition for P{a, a') reads 

PooPii = pPoiPio (A. 10) 

Note that Pier, a') is not normalized for a finite system and therefore a normalization con- 
stant J\f has to be taken into account. 



APPENDIX B: MAPPING TO AN EQUIVALENT ISING-MODEL 

Introducing Ising- variables cXi = ±1 instead of the lattice-gas variables Tj = 0, 1 (cij = 
2Ti — 1) one can look at the steady state as the equilibrium distribution P{cr) ~ e~^^^^^ 
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of an Ising model with Hamiltonian H — — J crjcrj+i + hY.iO'i In order to determine the 
couphng constant J and the external field h one interpretes the local probabilities P^.o-' as 
transfer matrices Pa,a' = ae~'^'^"''~'^^"'~^'^'^^'^, {P = !)• From eqn. (A. 10) it follows directly that 
e^"^ = p or J = log(p)/4 < 0. Therefore the corresponding Ising model has antiferromagnetic 
interactions. According to eqns. (A. 5) and (A. 6) one has in addition 



Dividing the two expressions on both sides of the two equations gives one eqn. without the 
constant a to determine the external field as 



Therefore the steady state corresponds to an Ising-model with antiferromagnetic (repulsive) 
interactions. Due to the conservation of the total number of particles (cars) one has to 
impose the constraint of a fixed magnetization to the Hamiltonian. The normalization is 
then simply the partition function calculated under this constraint. 



1 = 2 • Poi + + Poo 




p = Poi + Pii- 




^l-4(l-p)p(l-p)-(l-2p) 
2v^(l - P) 



(B.3) 
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TABLES 





s.bit (F77) 


particle (F77) intermed. (C) netw. (C) 


Sparc 10 


4.0 MUPS 
30 000 km 


3.4 MUPS 


1.9 MUPS 
14 000 km 


1.2 MUPS 
8 800 km 


PVM 
(5x SplO) 


19.0 MUPS 
140000 km 




8.9 MUPS 
65 000 km 




SX-3/ll(i) 533 MUPS 
1 CPN 4 000 000 km 




2.8 MUPS 

21 000 km 

_l www xv±X± 




GCel-3 102 MUPS 
1024 CPNs 750000 km 


211 MUPS 
1 550 000 km 


121 MUPS 
900 000 km 




iPSC 
32 CPNs 


83 MUPS 
630000 km 




35 MUPS 
260000 km 




Paragon 
64 CPNs 








290 000 km 


CM-5(i) 
32 CPNs 


173 MUPS 

1 300 300 km 




30 MUPS 

220000 km 




CM-5(i) 
1024 CPNs 








[> 1.7 e 6 km] 



CPN(s) has/have vector units (SIMD instruction set) 
using data parallel Fortran (CMF) 
(•^^ using message passing (CMMD) 
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TABLE I. Computational speed of different algorithms and different machines. The machines 
are a SUN SparclO Workstation, 5 coupled Workstations SparclO under PVM, a NEC SX-3/11 
traditional vectorcomputcr, a Parsytec GCel-3 massively parallel computer with 1024 relatively 
slow T805-CPU's, an Intel iPSC/860 Hypercube with 32 nodes, and Intel Paragon XP-IO/S with 
64 nodes, and a Thinking Machines CM-5 with 32 nodes, each node containing one Sparc processor 
and 4 vector units. — "s.-bit" refers to the single-bit coded site-oriented algorithm, "particle" and 
"intermed." to the particle-oriented and the intermediate one, and "netw." refers to a network 
implementation of the freeway network of Germany. — All values are valid for a vehicle density 
of c = 0.1 and system size of L = 1333 333 sites, corresponding to 10 000 single lane kilometers. 
MUPS corresponds to Mega- Updates Per Second, i.e., the number of sites updated per second 
divided by 10^. The other number is the extrapolated real time system size, i.e. the extrapolated 
system size (assuming linear speed-up) where the computation would be just as fast as reality. 
Values in brackets [] are estimates. 
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FIGURES 

FIG. 1. Evolution of different automata rules from an ordered initial state of density c = 0.5. 
Particles are moving to the right. Left: Same as (stochastic) asymmetric exclusion, except that the 
update is parallel. The same plot for exact asymmetric exclusion looks similar. Note that waves 
do not move in space. Right: Same as asymmetric exclusion, except that the maximum velocity 
is 5. Waves are now moving backwards, indicating a density above maximum flow. — Neither 
the parallel update nor the higher maximum velocity alone are sufficient to change the qualitative 
dynamics of the asymmetric exclusion model. 

FIG. 2. Evolution of our traffic model (maximum velocity Vmax = 5, parallel update) for 
density c = 0.3 (left) and c = 0.1 (right) from ordered initial conditions. Random sequential 
update with Vmax = 5 at density c = 0.3 (not shown) looks qualitatively similar to Fig. |l| left, 
whereas at density c = 0.1 and random sequential update the waves are moving to the right (not 
shown). The density of maximum throughput is much lower for the parallel update, and instead of 
the waves switching from backward to forward motion at this point, they tend to vanish completely 
(cf. ||2^). Moreover, the right picture illustrates the instability mechanism similar to the continuous 
model of [||. 

FIG. 3. (a) Fundamental diagrams flow j vs. density c for maximum velocities Vmax = li 2, . . ., 
5, and for a different fluctuation parameter p = 0.25 instead of 0.5 at Vmax = 5. (b) Fluctuations 
of local speed as a function of density {vmax = 5 and p = 0.5). 

FIG. 4. Partial densities co(c), ci(c), . . ., 05(0) (i.e. for velocities to 5) for mean field approx- 
imation, p = 0.5 

FIG. 5. Flow for Vmax = 00 in mean field approximation. From bottom to top, the randomiza- 
tion parameter p is 0.1, 0.3, 0.5, 0.7, and 0.9. 

FIG. 6. Convergence of the n-cluster approximations to the simulation result for the case 
Vmax = 1 and p = 0.5. Already the 2-cluster approximation is exact. 
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FIG. 7. Convergence of the n-cluster approximations to the simulation result for Vmax = 2 and 
p = 0.5. Already the 5-cluster approximation gives a good fit of the simulation data. 
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